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Abstract 

In this paper, we propose the SPR (sparse phase retrieval) method, which is a new phase retrieval 
method for coherent x-ray diffraction imaging (CXDI). Conventional phase retrieval methods effec- 
tively solve the problem for high signal-to-noise ratio measurements, but would not be sufficient for 
single biomolecular imaging which is expected to be realized with femto-second x-ray free electron 
laser pulses. The SPR method is based on the Bayesian statistics. It does not need to set the ob- 
ject boundary constraint that is required by the commonly used hybrid input-output (HIO) method, 
instead a prior distribution is defined with an exponential distribution and used for the estimation. 
Simulation results demonstrate that the proposed method reconstructs the electron density under a 
noisy condition even some central pixels are masked. 

1 Introduction 

X-ray free electron lasers (XFELs) can potentially provide us a novel mean to determine the three- 
dimensional (3D) structure of biomolecules from the diffraction images of single molecules [T,'2|. Con- 
ventionally, x-ray crystallography has been the principal tool to determine high-resolution 3D structures 
of proteins, nucleic acids, and their complexes. However, a critical process is the crystallization, where 
a sufhciently large crystal must be prepared. It is known that about 40% of biomolecules, particularly 
membrane proteins, do not crystallize. The coherent x-ray diffraction imaging (CXDI) by XFELs does 
not require a crystal and will change the situation. 

In order to realize the single molecule imaging, a straightforward scenario has been proposed [3111] 
and demonstrated with single mimivirus particles [S] . Firstly, a series of single molecules with the same 
structure are injected into vacuum and exposed to a very short (less than 5 fs) pulse x-ray beam. The 
diffraction image of the molecule with unknown orientation is recorded before radiation-induced structural 
degradation. This process is repeated until a sufficient number of images are recorded. The diffraction 
data are then classified into groups according to the orientations of the molecule and averaged to improve 
the signal-to- noise ratio [HIT]. Next, the relative orientations of the classes are determined to assemble 
a 3D coherent diffraction pattern in reciprocal space. Finally, the 3D electron density of the molecule is 
obtained by a phase retrieval method. 

The above scenario was proved to be effective for the 3D reconstruction of mimivirus particles [S]. 
However, the condition for a single molecule CXDI will be more severe since the number of the photons 
received by each pixel will be very limited. Some methods of reconstructing 3D diffraction intensity have 
been proposed for this problem [SHTU] . They are likely to be computationally intensive and the feasibility 
remains to be seen. 

Another scenario might be possible. After collecting a sufficient number of diffraction images with 
unknown orientations, 2D electron density images are first reconstructed from the 2D diffraction data, 
followed by the 3D reconstruction method used in transmission electron microscopy [11 . In cither sce- 
nario, the difficulty is due to the small number of photons scattered by a single molecule even though 
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x-ray beams that linac coherent Ught source (LCLS) has started to commission give us ihumination which 
is around 10® times brighter (^10^® photons/pulse/mm^ in the unfocused beam) than that of current 
synchrotrons [12] . 

We propose a new method for retrieving phases of a diffraction image which wiU be obtained by 
XFELs. The method is based on the maximum a posteriori (MAP) estimation of the Bayesian statistics. 
The posterior distribution is computed from the hkehhood and the prior. For the present problem, the 
Poisson distribution is appropriate for the hkehhood. The prior reflects the assumption of the electron 
density. Since the electron density is known to have a lot of zeros, we employ the exponential distribution 
which encourages the MAP estimate to have a lot of zeros (this is called a sparse solution [T31II1]. A 
similar idea is employed in [IS|). Therefore, we call our method, sparse phase retrieval (SPR) method. 

This combination is effective for the problem and is distinct from other Bayesian approaches proposed 
for different phase retrieval problems [inilll]. Compared to the widely used hybrid input-output (HIO) 
method p^H2Q] . the SPR method does not require the support region and is robust against the limited 
photon counts and the missing central pixels. In this paper, we focus on 2D phase retrieval problem. The 
resulting SPR method can be used directly for the second scenario. If 3D diffraction data are available, 
it is also possible to apply the SPR method. But reconstruction of 3D diffraction data from 2D images 
with unknown orientations is beyond the scope of the present paper. 



2 Proposed method 

2.1 Diffraction image and noise 

We first formulate the phase retrieval problem of a single diffraction image. Let f{x, y) > be the electron 
density of a molecule projected onto a two-dimensional (2D) plane. The coordinate is discretized, that 
is, y = 1, • • • , M . We introduce the following notation 

/ = {f{x,y)} - = {/n, • • • , /mm}, F = {F[u,v)} = = {i^u, • • • ,Fmm}, 

where, F is the Fourier transform of / defined as follows, 

Fu. = {F{f)U = ]g E f-y exp( ^ ■ (1) 

The coefficient of the Fourier transform differs from the commonly used 1/JVP. We prefer the above 
definition because the power of F and / is the same. We assume that the frame size is sufficiently larger 
than the molecule in the image and many components of / are zero. 

In the ideal situation, each pixel of a diffraction image receives the number of photons proportional to 
the power spectrum |F„^P, where {u,v) is the 2D index of the pixel. However, the XFEL measurement 
of a single molecule diffraction will be different because the flux of the x-ray laser is not sufficiently large. 
The number of the photons detected by each sensor will be small and behave stochastically. 

Let Nuv be the number of photons detected by the sensor at (u, v). We study the distribution of Nuv- 
The total number Naii = J2uv ^™ ^ stochastic variable which follows the Poisson distribution. The 
expected value of Naii is proportional to the intensity Ix of x-ray. Assuming each iV,j„ is independent 
and the distribution belongs to the same family, Nuv also follows a Poisson distribution. Let the expected 
value of Nuv be Suv and the distribution of Nuv becomes 

/AT ir. ^ Suv" — Suv) 
P{Nuv\Suvj = 7J—. • 

We only consider this counting process and do not consider other observation noise. Approximately, Suv 
is denoted cos'^ 6*, where a is a positive constant and 9 is the scattering angle (see Fig.[T]). 

Because is a function of u and v, we further rewrite it as a\Fuv\'^Cuv, where Cuv — c{u, v) = cos'^ 9. 
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Figure 1: A schematic picture of tfie single biomolecular diffraction imaging witli XFEL. 



The Fourier transform, F , is a linear function of /, and the probability of AT, where N — {Nuv\ — 
{iVii, • • • , Nmm], is denoted as follows 

uv 

Note that the scaling of l-Fu^l^ cannot be determined only from the observed scattering pattern and a is 
set to 1. This scaling indeterminacy is ignored here. It would be recovered from other knowledge, such 
as the total energy. 

In the following, we consider the density reconstruction as an estimation of /. Commonly used 
estimate is the maximum likelihood estimate (MLE) which maximizes p{N\f). The likelihood in Eq.® 
is maximized by setting = Nuv/cuv, but this is not sufficient for the estimation of / because the 

phase is lost. Some additional information is needed. 

2.2 SPR method 

Instead of the MLE, the MAP estimate of the Bayesian statistics is used for the SPR method. The MAP 
estimate maximizes the posterior distribution, where the posterior distribution is defined as the product 
of the likelihood function and the prior distribution. Similar framework has been proposed for different 
phase retrieval problems [inilTT] . 

The likelihood function is defined as a Poisson distribution as in Eq.([2]) and we discuss the prior 
distribution in this subsection. For a 2D image of a single molecule diffraction, it is natural to assume 
that the electron density is high in the central part and becomes zero in the peripheral part. But we do 
not know which pixel is zero. The prior should reflect this knowledge. 

It is known that the prior distribution with an exponential distribution encourages the MAP estimate 
to have a lot of zeros [T31[T3] and we use this for the SPR method 

Pif) = Il^^(-^^y)' Pif^^v) = Pxyexp{-pa:yfxy), fxy > 0, Pxy > 0, 

xy 

where each f^y is assumed to be independent. The hyperparameter p^y reflects the prior belief of fxy — 0. 
Larger the value of p^y, stronger the belief. In order to reflect the prior knowledge, we define the following 
one parameter function p{p)xy 

\( l + M\-^ ( 1 + M\21 
PKiAxy = pwxy, Wxy ^ai[x — J +[y — J >+b. (3) 
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Figure 2: The 3D plot of Wxy This reflects the prior knowledge of the electron density. 



The two parameters a and b were adjusted to make Wxy equal to at the center and 1 at the corners 
(see Fig. [2]) and /i adjusts the overall strength of the prior. We have employed a parabolic function, but 
if there is any prior knowledge, p{ij)xy should reflect it. 

The posterior distribution of / observing N has the following relation, 

p{f\N)<xp{N\f)p{f). 

It is convenient to take the logarithm of p{f\N) in order to compute the MAP estimate. By taking the 
logarithm and collecting terms related to /, the following function is obtained 

e^{f\N) = ^(Ar™ln|F„„|2 - jF^^pc™) - ^ ^ w.j,/,^,. (4) 

uv xy 

The first summation corresponds to the likelihood with the Poisson distribution (similar idea is used 
in |16)). the second summation represents the prior with the exponential distribution, and /i controls the 
balance between two terms. This combination is the key idea of the proposed method. The SPR method 
is summarized as follows. 

maximize £^{f\N) subject to fxy>0 {^<x,y<M). (5) 

We start with a large /i and decrease it. The MAP estimate is computed for each /i, and the best value 
of fji is selected according to some criterion. This issue will be discussed later. 

We note that our formulation is related to the compressed sensing (CS) [21], and the idea of CS 
has been applied to the phase retrieval problem [15]. Our approach is different from the work since the 
Poisson likelihood is included. 



2.3 Algorithm 

A gradient based algorithm is used to compute the MAP estimate. The algorithm is shown below starting 
with the following relation, 



— Ref F„^,exp(^ 



By defining the inverse Fourier transform as 

2'iri{ux + vy) 



(^-i(F)),, = ^^^^„„exp(- 



M 

UV 

the derivative of i^{ f\N) becomes 

^^il^^2Re{T-\g{F;N)))^^^Wxy, 

CJxy " 
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where g{F; N) = {guv{Fuv; Nuv)} and 

a (F -N ) - ( c \f 

\\-fuv\ J 

In order to compute the derivative of £^{f\N), T is apphed to /, then T^^ is apphed to g[F\ N). Thus, 
the computational cost is similar to the HIO method [19 . 

We employed a naive iterative updating rule for the MAP estimate of fxy, x,y = 1, ■ ■ ■ , M. 

/*+i^max(0, fiy + m ^^i^.T^ ). (6) 

where t is the iteration number. The "max" operation keeps fxy nonnegative. The positive variable r]t 
controls the step size and a simple line search of rjt accelerated the convergence largely. 

The initial value is generally important for gradient based algorithms. But for the SPR method, the 
influence is effectively reduced by controlling /i. The MAP estimate is first computed for a large /i and 
then for smaller /i's successively. For a very large /i, most components of the optimal / are and the 
influence of the initial value is negligible. The MAP estimate is then used for the initial value of a slightly 
smaller fj,, which is a good starting point in general. 



2.4 HIO method and Bayesian framework 

In this subsection, we discuss the relation between the widely used HIO method and the SPR method. 

The HIO method 18-20^ exploits the knowledge that the positive density is localized in a small 
support region. Let us deflne the support region 7, and the method is written as follows, 

minimize I c^/^ _ N^^^ ^ subject to fxy > 0, {x,y) £ 7 and fxy = 0, {x,y) i 7. (7) 

uv 

This problem can be formulated as the MAP estimation problem by setting the Gaussian likelihood 
function and the uniform prior distribution, that is, 

^ ' uv ^ 

where C is a positive number. If fxy is unbounded, above Puifxy) is improper. The MAP estimate is 
computed by solving Eq.©. 

Formally, the difference between the HIO method and the SPR method is the likelihood and the 
prior. But more importantly, the proposed SPR method does not require the support region. In the SPR 
method, a lot of components of the estimated electron density / become automatically. 



3 Numerical Experiment 
3.1 Simulated data 

The SPR method was tested with simulated diffraction data. The electron density of lysozyme, a 14 kD 
protein, was projected onto a 2D plane. The projected density was then Fourier transformed into recip- 
rocal space (m, v) to simulate a diffraction image. The density and the ideal diffraction image are shown 
in Figs. andlSb, respectively. 

Experimental diffraction images were simulated by converting the diffraction intensity into the number 
of photons per effective pixel {\/ aVf' according to the Poisson distribution with the mean, Suv, defined 
as follows. 
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(a) 2D electron density of a lysozyme protein. 







(b) A noiseless diffraction image 
with x-ray fluxes / = 5.0 X 10^^. 
The scale of the image is propor- 
tional to log(5„„ + 1). 




(c) Simulated diffraction image 
with x-ray fluxes / = 5.0 X lO^-"^. 
The total number of photons is 
4630. The scale of the image is pro- 
portional to log(iVu„ -|- 1). 



(d) Simulated diffraction image 
with x-ray fluxes / = 1.0 X 10'^-'^. 
The total number of photons is 957. 
The scale of the image is propor- 
tional to log(Af„„ + 1). 



Figure 3: Simulated 2D electron density (a), the diffraction image without noise (b), and simulated diffrac- 
tion patterns (c,d) of a lysozyme protein are shown. Each image is 308x308 pixels. The spatial frequency 
of the diffraction patterns is 0.0812 nm~^ and the unit of the x-ray fluxes I is photons/pulse/mm^. 
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Here, / is the incident x-ray flux (photons/pulse/mm^), Vc is the classical electron radius (2.82x10^^^ mm), 
A is the x-ray wave length (0.1 nm), cr is the oversampling ratio per one dimension (2), and L is the 
molecule diameter (6.16 nm) (see Fig. [T]). The Ewald sphere curvature was ignored but this is a good 
approximation for small angles. 

Two examples for scattered photons of x-ray fluxes with / = 5.0x10^^ and 1.0 x 10^^ photons/pulse/mm^ 
are shown in Figs. I^t and |3Jl, respectively. Although the fluxes in the above simulation are around 100 
times larger than that of the current LCLS, the sizes of the molecules of interest are generally larger than 
lysozyme. When the size is around 1000 kD, the total number of photons arriving at the 2D plane will 
be in the similar range as Figs. [3t and[5Ji. 

3.2 Results 




(a) Reconstructed electron density (b) Reconstructed electron density (c) Reconstructed electron density 
with fi = 10000. with fi = 100. with fi = l. 

Figure 4: Reconstruction of the electron density from the ideal diffraction image by the SPR method 
with different values of /i. 

The SPR method was applied to the data to reconstruct the density. Figure |4] shows the results of the 
SPR method applied to the ideal diffraction image. As ^ decreases, more components become positive. 
The reconstructed electron density becomes similar to the true density when ^ — 1. Thus, it is important 
to select the best /i. A simple strategy is to rank the results with a criterion. Many criteria have been 
proposed in statistics, such as, Akaike Information Criterion (AIC) and Bayesian Information Criterion 
(BIC) [22], but here we use a different one. Let /(/i) be the maximizer of £fj_{f\N) and the following 
error is used as the criterion 

Error^^ (/^) = ^nvilFif^Ulcuv'^^ - N^.'^^)^ ^ 

where F(fi)uv = J'ififJ-))uv This error function has been widely used for the phase retrieval problem and 
is employed here as well in order to rank the results. Roughly speaking, the error decreases as fj, shrinks 
since it controls the balance between the prior and the likelihood. But the error does not necessarily 
decrease monotonically because the log posterior in Eq. (j4]) is different from the error. Thus, we choose 
the fi which minimizes Errors? (/x). 

Figure [5] shows the Errors? (/i) for diffraction images in Fig. [31 The value of /i was varied from 1x10'' 
to 1x10^^. Note that although the horizontal axis of Fig. [5] increases from left to right, the value of ^ 
was decreased in the experiments. The best /i was selected for each diffraction image by Error i?(/i) and 
the reconstructed electron densities are shown in Figs. [6^, [Sb, and [St. 

We also applied the HIO method for the data. The central square (154 x 154 pixels) was set as the 
support region. The method has a problem with convergence, and we stopped the iteration loop when the 
incremental difference of the estimated density becomes small enough. Figure [7| shows the reconstructions 
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Figure 5: The values of Errori?(/i) for simulated diffraction images in Fig.[3l The black circle and triangles 
on the curves correspond to the reconstructed density images of Figs. ID 




(a) Reconstruction from diffraction (b) Reconstruction from diffraction (c) Reconstruction from diffraction 

image Fig. 2b by tfie SPR method image Fig. 2c by tfie SPR method image Fig. 2d by the SPR method 

with fj, = 0.001. ErrorF is with fj, = 0.1. Errors is 0.209. with fj, = 0.05. Errors is 0.277. 
1.60x10-3. 



Figure 6: Reconstruction of the 2D electron density by the SPR method. 




(a) Reconstruction from diffraction (b) Reconstruction from diffraction (c) Reconstruction from diffraction 
image Fig. 2b by the HIO method. image Fig. 2c by the HIO method. image Fig. 2d by the HIO method. 
ErrorF is 3.30x10-3. Errors is 0.220. Errors is 0.291. 



Figure 7: Reconstruction of the 2D electron density by the HIO method. 
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Figure 8: The values of Errori?(/i) for simulated diffraction images in Fig. |3] without the central pixel. 
The black circle and triangles on the curves correspond to the reconstructed density images of Fig. |9l 

obtained by the HIO method for the diffraction images in Fig. [3] The reconstructed electron densities 
are similar to Fig. [51 but the SPR method obtained better Errori?(/i) for all of the diffraction images. 

We note the computational time of the SPR method. The method was implemented with C and run 
on a desktop machine with an Intel Xeon X5570 (2.93 GHz 4 cores). The optimal electron density was 
computed for all the values of /x's (33 different values of fj,) within 10 minutes. 



3.3 Missing centers 

One of the problems of the phase retrieval from diffraction images is missing centers. The intensities 
around the center cannot be measured experimentally because the scattered beam is superimposed by 
the direct beam around the center and a beam stop is used to block the direct beam. 
The SPR method can be applied to the case easily by modifying eq. ^ as follows, 

uv^a xy 

where a is the active pixel set where the central pixels are excluded. The error function is also modified 
as 

Errori.(^) = =^ . (10) 

Firstly, we masked the central pixel from each diffraction image in Fig.[3]and applied the SPR method. 
By masking the central pixel, Fig.lJb lost around 5.7% of the photons, Fig.[3J; lost 294 out of 4630 photons 
and Fig. |3Ji lost 56 out of 957 photons. Figure [S] shows the errors. Note that the error values are better 
than those in Fig. [Sj but it does not mean the reconstructed densities in Fig. [HI are better than those in 
Fig. [6l because the errors are computed only for the active pixel set. Figure [TOl shows the results of the 
HIO method. Since the central pixel corresponds to a scaling factor of the electron density of the protein, 
the reconstructed images are quite similar to those reconstructed from the complete diffraction images. 
We see that the error values of the SPR method are better than those of the HIO method. 

Next, we further masked the central 3x3 pixels from each diffraction image and applied the SPR 
method. Note that Fig. lost around 35% of the photons. Fig. [5J; lost 1595 out of 4630 photons and 
Fig. [3ji lost 312 out of 957 photons. More than 30% of the photons are lost and the results are expected 
to be degraded largely. Figure [11] shows the errors. The reconstructed densities are shown in Fig. [T^l 
Figure [T3| shows the results of HIO method. We see that the error values of the SPR method are better 
and the reconstructed results are more stable than the HIO results. It is worth noting that the SPR 
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(a) Reconstruction from diffraction 
image Fig. 2b without the central 
pixel by the SPR method with /i = 
0.0001. Errorj? is 1.78x10"^. 



(b) Reconstruction from diffraction 
image Fig. 2c without the central 
pixel by the SPR method with /i = 
0.002. Error is 0.156. 



(c) Reconstruction from diffraction 
image Fig. 2d without the central 
pixel by the SPR method with /i = 
0.005. Error j7 is 0.195. 



Figure 9: Reconstruction of the 2D electron density by the SPR method. 




Figure 10: Reconstruction of the 2D electron density by the HIO method. 



method worked well even if the central 5x5 pixels, a half area of the centrospeckle pattern, were masked. 
In that condition, 54.4 % of the total photons are lost. 



4 Discussion 

We have proposed a new phase retrieval method. The SPR method is based on the MAP estimation of 
the Bayesian statistics, it will be effective for single molecule diffraction images which will be obtained 
by XFELs. 

From a Bayesian viewpoint, the difference between other phase retrieval methods including the HIO 
method and the SPR method is summarized in the likelihood and the prior. Many methods use the 
squared loss, which corresponds to a Gaussian likelihood, while the present method employs a Poisson 
distribution. This is suitable for CXDI data because the number of the scattered photons detected at 
each pixel is a counting process. Similar idea has been found in a related work [16] . The prior in HIO 
method corresponds to a uniform distribution with a bounded support, while the proposed method uses 
the exponential distribution. Similar prior has been widely used in statistics [131 . This combination 
shows a new promising direction for phase retrieval. 
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Figure 11: The values of Errori?(/i) for simulated diffraction images in Fig. |3] without the central 3x3 
pixels. The black circle and triangles on the curves correspond to the reconstructed density images of 

Figs, m 




(a) Reconstruction from diffraction (b) Reconstruction from diffraction (c) Reconstruction from diffraction 

image Fig. 2b without the central image Fig. 2c without the central image Fig. 2d without the central 

3x3 pixels by the SPR method 3x3 pixels by the SPR method 3x3 pixels by the SPR method 

with II = 0.0001. Errovp is with = 0.002. Errorj? is 0.222. with = 0.005. Error is 0.267. 
2.34x10"^. 



Figure 12: Reconstruction of the 2D electron density by the SPR method. 




(a) Reconstruction from diffraction 
image Fig. 2b without the central 
3x3 pixels by the HIO method. 
Errors- is 0.128. 




(b) Reconstruction from diffraction 
image Fig. 2c without the central 
3x3 pixels by the IflO method. 
ErroTp is 0.487. 



(c) Reconstruction from diffraction 
image Fig. 2d without the central 
3x3 pixels by the HIO method. 
Error p is 0.557. 



Figure 13: Reconstruction of the 2D electron density by the HIO method. 
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The proposed method has been tested with simulated data. The electron densities are reconstructed 
with a reasonable computational cost, and errors of the results are smaller than the HIO method. We 
further tested the method with the diffraction data where some central pixels are masked. The proposed 
method is applied without any major modification, and the results are quite stable even some pixels are 
masked. 

The SPR method reconstructs the density from an incomplete observation utilizing the sparsity. This 
idea is similar to CS |21) . In |15] . CS is applied to the phase retrieval problem, however, in our problem, 
the number counts are very limited, and we needed to extend the model. 

Finally, we list some of our future works to improve the SPR method. One is the hyperparameter. 
The hyperparameter p{n)xy can be modified to reflect the knowledge of the true density. This may result 
in a better estimate. Another issue is the algorithm. The densities for the whole sequence of in Fig. [3] 
are computed with a simple line search algorithm but advanced optimization methods may speed it up. 
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